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Abstract 

Analytical derivations of stress intensity factors (SIF’s) for a multicracked plate can be com- 
plex and tedious. Recent advances, however, in intelligent application of symbolic computation 
can overcome these difficulties and provide the means to rigorously and efficiently analyze this 
class of problems. Here, the symbolic algorithm required to implement the methodology de- 
scribed in Part I is presented. The special problem-oriented symbolic functions to derive the 
fundamental kernels are described, and the associated automatically generated FORTRAN 
subroutines are given. As a result, a symbolic/FORTRAN package named SYMFRAC, capa- 
ble of providing accurate SIF’s at each crack tip, has been developed and validated. 

Simple illustrative examples using SYMFRAC show the potential of the present approach 
for predicting the macrocrack propagation path due to existing microcracks in the vicinity 
of a macrocrack tip, when the influence of the microcracks’ location, orientation, size, and 
interaction are accounted for. 

Nomenclature 

direction cosines between two local coordinate systems 
strain tensor 

offset of notch-microcracks system with respect to Y axis 
four roots of the characteristic equation 
Poisson’s ratio 

far-field and total stress field, respectively 
components of stress in global coordinate system 
angle defining orientation of local coordinate system 
nprmalized real variables 

Fourier transform of Airy stress function with respect to x variable 
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half crack length 

coefficients of strain-stress relationship (compliant matrix) 
normalized radial (tip) distance 
auxiliary functions 

mode-I and mode-II stress intensity factors 
Fredholm kernels 
normal traction at crack surface 
shear traction at crack surface 
Fourier variable 

position vector defining the origin of a local coordinate system 
components of the position vector rj 

displacement associated with the x and y coordinates, respectively 
weight function 
local and global coordinates 
kernel matrix 

functions of s in Fourier space (i.e., constants in x, y-real space) 

Young’s modulus 
Airy stress function 
discrete auxiliary function 
loading vector 

1 Introduction 

The computer has become an indispensable tool for both engineering analysis and design. 
Advanced computing techniques provide powerful tools for computationally-intense applications in 
engineering. Symbolic computation specializes in the exact computation with numbers, formulas, 
vectors, matrices, equations, and the like. Numerical computation, on the other hand, uses floating- 
point numbers, and approximates computations to solve problems of practical interest. The two 
approaches are complementary and, when combined into an integrated form, can be very powerful 
for engineering applications. 

Analyzing the interaction of microcracks analytically, as discussed in Part I (Binienda et al. 
1992), involves extensive manipulation of complex mathematical expressions. To date, the method- 
ology has been developed for analyzing multiple cracks within an isotropic material. A similar 
methodology can be implemented for anisotropic and/or nonhomogeneous materials, and for com- 
plex, nonstraight multicracked configurations. However, because of the complexities involved, it 
is impossible for researchers to investigate this general problem without relying on the power of 
symbolic computations. 

Presented here is a symbolic manipulation package (SYMFRAC, from SYMbolic FRACture) 
capable of providing both analytic derivation and FORTRAN code generation for n straight, fully 
interacting cracks in an isotropic plate that is subjected to in-plane loading. The immediate benefits 
that can be realized are (1) reduced tedium of manual derivation, (2) increased reliability of the 
derived equations and, hence, the final analysis results, and (3) improved numerical efficiency and 
accuracy for multiple interacting cracks. 

This paper begins with a description of the associated symbolic algorithm and is followed by 
subsequent sections describing in detail the three key steps involved in the underlying formulation 
given in Part I (Binienda et al. 1992). It concludes with two numerical examples that illustrate 
the major capabilities of SYMFRAC: (1) a four-crack problem in which three microcracks interact 
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with a larger notch and (2) a fully interacting multicrack problem, with two notches and eight 
microcracks. 


2 Computational Algorithm 

The general theoretical formulation for a multicrack mixed boundary value problem has been 
discussed in Part I (Binienda et al. 1992) and can be classified into four main steps: (1) the 
derivation of the local stress equations for each crack in its respective local coordinate system, (2) 
formulation of the total perturbation stress field for each crack, (3) formulation of the singular 
integral system of equations, and (4) the solution for the stress intensity factors of this singular 
integral system via the discrete auxiliary functions. 

Each of these basic steps involve numerous, tedious intermediate steps; this suggests that sym- 
bolic computations may be an attractive tool for automating the derivation and solution of this 
class of mixed boundary value problems. 

The required algorithm to accomplish such automation can be divided into the following 11 
steps. 


(1) Convert the governing equation, 


£i F{x ’ y)+ 




( 1 ) 


into an ODE by Fourier transform; that is, 


s 4 $ - s 2 <I>" + = 0 


(2) 


(2) Solve the ODE. For the isotropic case, the roots of the characteristic equation are two identical 
real roots; thus the general solution is 

{C[ + C' 2 y)t ya + (C' + C' 4 y)e- y ’ (3) 

(3) Use the inverse Fourier transformation and the condition that F(x,y) must go to zero at 
infinity, to obtain the following Airy’s stress function: 

F(x, „+) = + C^e-U'e-”) for y > 0 (4) 

Z7T 

and 

F(x, y~) = ^-((C 3 + C 4 y)el-l'e-“] for y < 0 (5) 

where Cj(j = 1, ...,4) are arbitrary functions of s. 

(4) Find the stress, strains, and displacements through differentiating, applying a constitutive 
relationship and integrating with respect to strain. 

(5) Introduce auxiliary functions /i(x) and fi{x) such that 

/i(*) = ^[u + (*,0)-u - (x,0)]. (6) 

and 

h{x) = — [u + (x,0) - iT(x,0)]. (7) 
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(6) Solve for C\,Ci, C 3 , and C4 in terms of the auxiliary functions /i(x) and / 2 (x) by using the 
continuity conditions for a yy and a xy ; that is, 

ff w( x >°) = ( 8 ) 

<r£,( x >0) = y( x >°) ( 9 ) 

(7) Find the final form by substituting the local stresses in terms of the auxiliary functions fi(x) 
and / 2 (x), and integrate with respect to the Fourier variable. 

(8) Find the total stresses 

n—1 

j<rJz( x j,yj) = *?*(*;> yj) + 5Z °T P z[ x i>{ x i, w). y P ( x ^y 3 )} ( 10 ) 

p=i 

by using coordinate transformations, 

r jX + xj cos <pj - yj sin <pj = r pX + x p cos y p - y p sin <p p (11) 

r jY + Xj sin ipj + yj cos tpj = r pY + x p sin v? p 4- y p cos <? p (12) 

and the stress transformations 

<7 tz = PlrPmzVlm (13) 

where /? (r ,/3 mz are the direction cosines for the (xfc,y*) and (x p ,y p ) coordinate systems. 

(9) Identify the Fredholm kernels (ker p ) and normalize them, where p = 1,2, 3, 4. 

(10) Apply a collocation technique and generate the FORTRAN subroutine for the discrete auxil- 
iary function vector {G}: 

{G} = (14) 

where [A] is a fully populated matrix of the coefficients obtained from the Fredholm kernels, 
and {7£} is the applied loading function. 

(11) Evaluate the stress intensity factors (SIF’s). 


2.1 Symbolic Algorithm for Local Stress Formulation 

The general form of the Airy stress function F(x, y) for a plane problem is given in terms of the 
four unknown coefficients Ci to C4 (see Eqs. (4) and (5)). Two of the four are dependent because 
of the continuity conditions. The remaining two may be expressed conveniently in terms of the 
auxiliary functions as shown in Part I (Binienda et al. 1992). Substituting Ci to C4, as described 
in Part I into the stress equation for y > 0 or y < 0, the stress equations become 
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In Eqs. (15) to (17) the range of integration is -a < t < a for the coordinate variable and 
-oo < s < oo for the Fourier variable. Integration with respect to the Fourier variable s can be 
separated into a negative portion (-oo,0] and a positive portion [0,+oo). Then the variable of 
integration can be substituted, with 5 ' = -5, and the order of integration switched such that all 
integration is over the domain [0, + 00 ). Hence, the stress equations are 
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It is important to note that at this stage the function integrate in MACSYMA (MATHLAB 
GROUP 1984) cannot be applied directly to perform the previous integration because of the oc- 
currence of an infinite loop. Instead, intermediate variables are introduced into the integrals and 
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the following recursive forms are applied in a heuristic manner to actually perform the required 
integration; that is, 

f°° 

h = / t pa ds where p = p(x,-y,t ) (21) 

J 0 


and 


d 

■n — 1 

dy 


In = ^In-l 


for n = 2,3, ... 


such that the result of the preceding integration (see Gradshteyn and Ryzhik 1980) is, 

(- 1 )" 


In = 


For the isotropic case of n=l and 2: 


P n 


for n = 1,2,3, ... 


h = 


1 


h = 


Therefore, the resulting stresses become 
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[: y 2 + (t - x) 2 ] 2 y 2 + (t- x) 2 
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y 2 - (t - x) 2 2 

[y 2 + (t-x) 2 ] 2 y 2 + (t-x) 2 
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The above evaluation of the Fourier integral is accomplished by invoking a single function called 
ISTRESS, thus making the details of this section transparent to the user. 
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2.2 Symbolic Algorithm for Fredholm Kernels 

As we showed in Part I (Binienda et al. 1992), prior to obtaining the Fredholm kernels, we 
must compute the total stresses due to all cracks within the plate (see step (8), Eqs. (10) to (13)). 
We accomplish this through both coordinate and stress transformations of all surrounding cracks 
to a chosen ] th crack as described by Eq. (10). This transformation and summation procedure is 
continued until the total stress state (i.e., G<rJ,) and G<rJ y )) at each crack surface has been found. 
Given the total stresses, we obtain the final form of the normalized Fredholm kernels by performing 
a variable transformation, (i.e., x=a£; t=ar; for -1< £ and r <1) and applying crack boundary 
conditions (i.e., y = 0). 

Hence, for the case of n cracks the whole system of singular integral equations is 

A = 4^7 O’-i keri ^ 11<fr + /-I ker 2/i2 rfr (29) 

+... + fli keri/( n _i)i<fr + ker 2 /(n-i)jdr + £ /ij ^drj 

Av = * {/-i kera/ndr + fl , ker 4 f 12 dr (30) 

+ ••• + fh ker 3 /( n _i)idr + fi j ker 4 f( n -i)idr 4- - f_ x dr | 

where the four distinct kernels (i.e., ker p for p= 1 ,2,3,4) are shown in the Appendix of Part I 
(Binienda et al. 1992). These four kernels can be translated into FORTRAN code directly through 
the use of a built-in command in MACSYMA. The resulting generated code is given in Appendix 
A at the end of this paper. Note that these kernel functions are already modified by the associated 
Lobatto-Chebyshev parameters discussed in the next section. 

2.3 Solution of Discrete Auxiliary Functions via the Collocation Tech- 
nique 

The Lobatto-Chebyshev collocation integration technique is used to transform the preceding 
system of singular integral equations (represented by Eqs. (29) and (30)) into a system of alge- 
braic equations. These equations combined with the single valued conditions, described in Part I 
(Binienda et al. 1992), result in the following 2nm system of algebraic equations: 

T[ Wjn{T ~ + ker,/ai(r r H + ker 2 / 22 (r r H 

r =i »(*V - 6) 

• • • + keri/ nl (r r )u) r + ker 2 / n2 (r r )ui r ] = 4a 1 i(i<rJ y ) 

for z = 1, ..., m — 1 (31) 

VI ™’W r z L + ker 3 / 2 i(T r )uv + ker,/«(T r )uv 

fZi*( T ' -6) 

• • • + ker 3 /„i (r r )ui r + ker 4 /n 2 (rr)u? r ] = 4a n (i<rJ y ) 

for z = 1, ..., m — 1 (32) 


A,. ... . , , / x , «v/2i(ty) 

2_j[ker i/ 2 i(r r )i^ r + ker 2 f 22 {T T )w r + 

• • + ker i/ n i(r r )ui r + ker 2 f n2 {r T )w T } = 4a„( 2 <rJ y ) 
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for z = — 1 


(33) 


X^ ker 3/2i(Ty)u>r + ker 4 / 22 (r r )u; r + 

r=l *{T r — 

• • • + ker 3 f nl (r r )w r + ker 4 / n2 (r r )u; r ] = 4a n ( 2 <7^) 

for 2 = 1, m — 1 (34) 


53[keri/ n (r r )uj r + ker 2 / 12 (r r )u; r + keri/ 21 (T r )iu r 

^r/nl(7- r ) 


r= 1 


+ ker 2 / 22 (r r )uv + 1- 


v { T r ~ 6) 

for z = 1 , m — 1 


] = 4a n ( n <7i.) 


(35) 


[ker 3 /n (r r )uj r + ker 4 / 22 (r r )tt> r + ker 3 / 21 (r r )iu r 

r=l 

+ ker 4 / 22 (r r )uv + h y-'— ] = 4an( n <rJ ) 

7r(r r -6) 


/or 2 = 1, m — 1 (36) 

m 

H/nKH = 0 (37) 

r = l 
m 

f\2{T r )w r = 0 (38) 

r = l 
m 

X)/2l( T r) U, >- = 0 (39) 

r=l 

m 

2 /22(r P )wr = 0 (40) 

P=1 


m 

= 0 (41) 

r= 1 
m 

X)/"2( T r)^ r = 0 (42) 

r=l 


where m is the number of collocation points and r r , to r and are the associated abscissae, weights, 
and collocation points described in Part I (Binienda et al. 1992). 
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The above equations can be rewritten compactly by using matrix notation; that is, 

m{g> = m ' ( 43 ) 

where the kernel matrix [A] is shown in schematic form in Appendix B and the FORTRAN code to 
assemble this matrix is given in Appendix C. Similarly, the associated FORTRAN code to assemble 
the loading vector {7£} is also shown in Appendix C. 

Although, only two key FORTRAN subroutines have been shown here, an entire FORTRAN 
code named SYMFRAC has been developed. This code has the capability to calculate the stress 
intensity factors at each crack tip as defined in Part I (Binienda et al. 1992) for an isotropic 
plate subjected to an in-plane stress field with n cracks of arbitrary geometry. This code has been 
employed to give the numerical results in this paper and in Part I (Binienda et al. (1992)). 


3 Numerical Application 

In this section two problems of a cracked, brittle infinite isotropic plate are studied. The first 
problem is composed of a notch and three interacting microcracks. We will obtain mode-I and mode- 
II SIF’s for the inner crack tips under two loading cases: normal and shear far field stress states. 
The geometry of this problem lacks any symmetry and represents a physical problem involving the 
influence of a cloud of microcracks on the fracture zone of a major notch. The second problem 
represents a symmetric crack configuration composed of two interacting sets of a notch and cloud of 
microcracks similar to problem 1. Here, however, each notch is associated geometrically with four 
interacting microcracks. This problem, which is also physically possible, illustrates various cases of 
propagation of two major notches that are influenced by their mutual interaction and, in addition, 
by the existence and interaction of a cloud of microcracks in front of these notches. Symmetry with 
respect to the origin of the global coordinate system is assumed in order to simplify the graphical 
reporting of SIF’s obtained for all inner and some outer crack tips. 

Macroscopically, the extension of the notch can be simulated through the connection with a sur- 
rounding microcrack. The specific microcrack involved, however, depends on the loading conditions 
and crack interaction effects (which are dependent upon the particular location, size and orientation 
of the surrounding microcracks). After identifying this specific microcrack we can assume that the 
microscopic extension of all cracks will occur in a self-similar manner. Consequently, the fastest 
growing microcrack will connect first with the notch, to create a large kinked or branched macro- 
crack (notch). The criterion for microcrack propagation is assumed to be represented by a critical 
SIF value. For the sake of illustration, let us also assume that the total critical SIF is a constant 
material property such that the maximum total SIF obtained by using mode-I and mode-II of the 
SIF’s (kj and k 2 ) normalized with respect to the far field stress (<r? p ) and y/a^. Thus for the j th 
crack tip, 

*Li = + (^)’U» < 44 ) 

is used as a crack selection criterion for propagation. 

Clearly, alternative criteria for crack propagation can be selected; therefore to maximize the 
future utility of the present results, the SIF’s are reported separately. Note, that the objective 
of these examples is not to examine the validity of a particular criterion, but merely to illustrate 
qualitatively the capabilities of SYMFRAC. 
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3.1 Interaction of One Notch With Three Microcracks 

Consider the problem of the extension of a horizontal notch under two far field loading conditions; 
that is, a normal and a shear stress. Assume that within the fracture zone of the central notch, a 
cloud of microcracks has developed due to localized flaws, grain boundaries, and/ or other fabrication 
and material factors. For the sake of simplicity, consider only three microcracks that are situated 
radially with respect to a horizontal notch, as shown in Fig. 1. The plate is subjected to the 
previously mentioned two loading conditions: Case 1 — a normal far field stress, and Case 2 a 
shear far field stress. 

3.1.1 Case 1: Normal stress, <tyy 

The variations of ki and ik 2 for the inner crack tips as functions of the normalized distance 
d is shown in Figs. 2 and 3. The distances d \ 2 , d\ 3 , and dn measured between the inner tips 
of the notch and microcracks 2, 3, and 4, respectively, are taken to be equal to each other (t.e., 
d = di 2 — di 3 = d 14), and they are normalized with respect to half the notch length, a\. The 
associated crack lengths, normalized with respect to a 1( are ai = 1 and a 2 = <13 = <*4 = 0.1. 

Figures 2 and 3 illustrate the mutual influence of the notch and microcracks on the mode-I 
and mode-II SIF’s, respectively. Note that the significant influence of the notch on the SIF’s of all 
microcracks begins once d < 0.1; however the notch continues to control (in that it will propagate in 
a self-similar manner), provided that the initial radial (tip) distance d of all microcracks is greater 
then 0.03. If d < 0.03, the mode-I SIF for crack 3 (oriented at 45°) becomes the largest, (see 
Fig. 2) Similarly, mode-II SIF for crack 3 dominates for all distances d, but in general it is much 
smaller than mode-I SIF’s. The total SIF is the largest for crack 3 when d is approximately less 
then 0.075; therefore the 45° crack will grow to connect with the notch such that it will become 
a kinked macrocrack. It should also be noted that for this loading condition when d approaches 
infinity, mode-I SIF asymptotically becomes yj ajcos 2 # and mode-II SIF becomes yJajsinOcosO. 

Now consider the influence of the size of inclined microcracks (t.e., cracks 3 and 4). It can be 
observed from Figs. 4 and 5 that the length of microcrack 3 will accelerate self-similar crack propa- 
gation of the horizontal notch (due to a maximum ki) until a 3 becomes larger than 0.95, whereupon 
k total of crack 3 dominates. Thus, kinking will occur in the direction of crack 3. Alternatively, Figs. 
6 and 7 show that self-similar extension of the notch is independent of the size of crack 4, the 90° 
crack. In fact, Fig. 6 shows that crack 4 shields the notch by reducing its SIF’s once a 4 > 0.8. 

Note, that in Figs. (4) to (7) the centers of the microcracks are maintained at r 2 = r 3 = r 4 =1 
(where r } = dj + aj ), while the length of the particular crack under consideration varies from 0 to 
1 and the size of the other microcracks remain constant at aj = 0.1. 

3.1.2 Case 2: Shear stress, <?xy 

Next, a similar parametric study is conducted for the case when a far field shear stress is applied 
to the plate. As presented previously the variation of k\ and /: 2 with respect to d for the notch as 
well as all three microcracks is shown in Figs. 8 and 9, respectively. Clearly, the notch will grow 
in a self similar manner under mode-II for most values of d; however, the notch may kink in the 
direction of crack 3 if d is very small. The trend for mode-II, shown in Fig. 9, illustrates that for 
small values of d the SIF’s of cracks 2 and 4 are reduced due to shielding effects. The dominant 
mode of fracture for the notch is mode II when d > 0.03. 

Figures 10 and 11 show the mode-I and mode-II SIF’s for the inner crack tips when the size of 
crack 3 is increased, while the centers of all other microcracks are held constant at rj = 1. Results 
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indicate that for all lengths of crack 3, even those very large, the notch would not kink, but rather 
grow in a self-similar manner under mode-II conditions. Conversely, the SIF’s calculated when the 
size of crack 4 (i.e., 90°) increased ( see Figs. 12 and 13), indicate that the notch will grow (1) in 
a self-similar manner under mode-II conditions when a* < 0.63; (2) will kink in the direction of 
crack 4 when 0.63 < a 4 < 0.95 under mixed-mode conditions; or (3) will propagate once again in a 
self-similar manner, but now predominantly under mode-I, when a 4 > 0.95. 


3.2 Interaction of Two Notches With Eight Microcracks 

Consider the problem of two horizontal microcracks, possessing sharp notches, embedded in a 
plate parallel to the X-axis (see Fig. 14). The plate is subjected to a normal stress <ryy at infinity. 
Furthermore, assume that within the fracture zone of the notches, two clouds of microcracks have 
developed due to localized flaws, grain boundaries, and/or other fabrication and material factors. 
For the sake of simplicity, consider only four microcracks associated with each notch, such that 
they are situated radially with respect to each notch as shown in Fig. 14. The orientations of the 
microcracks within each cloud are 0°, 30°, 60°, and 90°, and the radial tip distance between the 
notch and its associated cloud of microcracks is d. 

The second parameter characterizing the geometry of the multiple crack system is A, which 
represents the offset of the system with respect to the global vertical axis Y. For the case of A > 0, 
the lower system of cracks (i.e., notch 1 and microcracks 3,4,5, and 6) is shifted to the right of 
the vertical axis, whereas the upper system (notch 2 and cracks 7,8,9, and 10) is shifted to the 
left. Note: The origin of the global coordinate system X and Y is always taken to be the point of 
symmetry of the crack configuration. In this way the presentation of the results can be simplified; 
however, this in no way implies any restriction, due to the symmetry, on the calculation of the 
SIF’s. Finally, the parameters d, A, and the half crack lengths aj ( j = 2,3,4,..., 10) are normalized 
with respect to the half notch length ai (where ai=a 2 ). 

For convenience, the location of each crack’s local coordinate system (i.e., center of the crack) 
is taken to depend upon the tip distance d, the offset A, and the half crack length dj. Hence, the 
centers of notches 1 and 2 are 

f\x = A — a.\ (45) 

r X Y = -{d + a 6 ) (46) 

and 

r 2 x = ~r\x 
r 2 y = —r x y 

and the centers of microcracks 3 through 6 and their counterpart microcracks (7 to 10) are 


rjx = A + (d + a,) cos ip } (49) 

rjy = r lY + (d + aj) sin <pj (50) 

r (13-j)X = ~ r jX (51) 

r (13 -j)Y ~ r jY (52) 


where rjx and rjy are the X and Y components, respectively, of the vector between the origin of 
the global and local coordinate systems and j =3, 4,. ..,10. 


(47) 

(48) 
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3.2.1 Influence of crack tip distance d 

Consider the case when oj s aj = 1, the offset A is equal to 0.1, and the half length of each 
microcrack aj = 0.1 (where j = 3,4,..., 10), and d varies from 0 to 0.5. The resulting inner SIF’s 
Jfcx and k 2 (representing mode I and mode II), are shown in Figs. 15 and 16, respectively. From 
examining the strength of k x and k 2 , it is clear that for all d values mode I is dominant and for 
d > 0.025 the inner notch tip would propagate in a self-similar manner. However, for d < 0.025, the 
Jki value for the 30° microcracks exceeds all other k x values, thus indicating propagation of the 30° 
microcracks and possible connection with the inner tip of the notch, thereby creating a macroscopic 
kinked crack. Although mode II (see Fig. 16) is significantly smaller in magnitude relative to mode 
I, and therefore does not play a role in the propagation of the various cracks, it is interesting to 
note that for d > 0.025 the 90° microcracks have a maximum k 2 , whereas for d < 0.025 the 60° 
microcracks do. Also, as expected when d approaches 0.5 (i.e., 1/4 of the notch length) the influence 
of the microcracks on the ki of the notch is minimal. Clearly there are two competing mechanisms: 
the first reduces k 2 because of crack shielding (observed at high d values), and the other increases k 4 
at small d values thereby simulating an increase in damage density at the notch tip that promotes 
crack propagation. 

3.2.2 Influence of offset A 

Here, as in the previous case, = a 2 = \ and aj = 0.1 for j = 3,4, ...10; but now d is held fixed 
at 0.1 while A is varied from -0.5 to +0.5. As before, the inner SIF’s (/►! and £ 2 ) a r e shown versus A in 
Figs. 17 and 18. Note, that all cases in which microcrack overlap occurs (i.e., —0.2 < A < 0.15) have 
been skipped. As in the previous section, two primary effects, one dealing with the “density” of the 
damage zone at the crack tip and the other with shielding are observed in Figs. 17 and 18. Clearly, 
in Fig. 17 as A approaches -0.2 from the left, the “density” of the damage zone (characterized by the 
number of microcracks per unit length in front of the notches) increases, consequently increasing 
the strength of k x . The second effect, again referring to k^ in Fig. 17, is attributed to the shielding 
for A > 0.15 in that each notch shields the other as well as their associated cloud of microcracks. 
Again, although mode II values are smaller than those of mode I, it is interesting to note that k 2 
of the notch becomes maximum when shielding occurs (i.e., 0.15 < A < 0.5). 

3.2.3 Influence of microcrack lengths a 4 and a 5 

Finally, the lengths of inclined cracks (e.g., 30° and 60°) are examined. Here, s is chosen to avoid 
crossing of the inclined crack with either notch. For example, consider the case in Figs. 19 and 20 
where a x = a 2 = 1.0, A = 0.8, d = 0.1, and a 3 = a 5 = a 6 = 0.1 while a 4 varies. Similarly, in Figs. 
21 and 22 <*! = a 2 = 1.0, A = 1.0, d = 0.1, and a 3 = a 4 = a 6 = 0.1 as a 5 varies. In these figures 
(Figs. 19 to 22) in addition to all SIF’s associated with the inner crack tips, the SIF’s associated 
with the outer crack tips of the notches and the varying inclined microcracks are shown. 

Comparison of Figs. 19 and 20 shows that mode I once again dominates. Therefore, let us focus 
our attention on the mode-I SIF for the outer tip of the notch. It is apparent that fcj=l for values 
of a 4 < 0.2, thus indicating no crack interaction at this tip and propagation away from the damage 
zone. However, when a 4 > 0.2, k 2 for the outer crack tip dramatically drops off due to shielding by 
the 30° inclined microcrack associated with the opposite notch. Now focusing our attention on the 
ki associated with the inner tip of the 30° inclined microcrack, we see that for a 4 > 0.2 the present k 2 
is maximum, thus indicating propagation of the inclined microcrack towards the associated notch; 
hence the potential kinking' of the inner notch tip. Figure 19 illustrates the importance of the initial 
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damage configuration and its impact on the mechanism of crack propagation. 

Comparing Figs. 21 and 22, which are associated with the case where the two notches are 
exactly aligned as shown in the inserts, we observe for the first time the potential dominance of 
mode-II SIF’s. When a 5 < 1.2, both inner and outer notch tips are shielded because fci < 1.0. 
However, the inner notch tip SIF is higher than the outer tip SIF because of the interaction of the 
associated cloud of microcracks. Some additional shielding of the outer notch tip occurs when the 
length of the 60° microcrack is increased up to approximately 0.4. 


4 Summary and Conclusion 

The symbolic algorithm and key FORTRAN subroutines for constructing and calculating the 
SIF’s, with the singular integral equation technique, for a multicracked mixed boundary value prob- 
lem have been presented. This work has resulted in the development of special problem-oriented 
symbolic functions (e.g., ISTRESS) running under MACSYMA that simplify and automate the 
derivation process. Also MACSYMA’s automatic FORTRAN generation capability has been uti- 
lized to produce the Fredholm kernel subroutines that constitute the core of the resulting FORTRAN 
portion of SYMFRAC. 

The accuracy of the present technique was confirmed with available solutions in the literature 
for a two interacting crack problem shown in Part I. Also, two fully interacting multicrack problems 
were studied here to illustrate the capabilities of SYMFRAC. 

Numerous observations where discussed using a simple crack propagation criterion. For example, 
(1) a notch-like crack can change its propagation direction (kink or branch) through connection with 
pre-existing microcracks; (2) the notch may propagate either toward or away from the surrounding 
cloud of microcracks; (3) two potential competing effects exist — shielding, which reduces the SIF s, 
and damage density, which increases the SIF’s. In all of the preceding examples, results depend on 
the size, orientation, and distribution of the interacting cracks. Thus, the most important conclusion 
demonstrated is that the current damage configuration and the loading history dictate the future 
damage growth. 

The power and usefulness of SYMFRAC is apparent in that the size, orientation and distribution 
of n fully interacting cracks in an isotropic plate can be studied. Furthermore, the methodology has 
now been established so that extension to more complex problems, in which anisotropic materials, 
kinked and branch cracks, and/or non-homogeneous materials are addressed, is straight-forward. 
Finally, it is apparent that the rigorous symbolic development of the system of singular integral 
equations and the associated automatic FORTRAN implementation is responsible for the speed 
and accuracy of the numerical calculations. 
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Appendix A - FORTRAN Code for the Discrete Fredholm Kernels 

**************************** **************** ******** 

* 

* Calculate krl kr2 kr3 and kr4 

* 

**************************************************** 

do 320 m=l,4 
do 310 k=l,ncr 
do 311 1=1, na- 
if (k .NE. 1) then 
r2x=rx(k) 
r2y=ry(k) 
gf2=gf(k) 
a2=al(k) 
rlx=rx(l) 
rly=ry(l) 
gfl=gf(l) 
al=al(l) 
gh=gf2-gfl 

pi =(r2y-rl y )*dsin(gfl )+(r2x-rlx)*dcos(gf 1 ) 
p2=(r2y-rly)*dcos(gfl)-(r2x-rlx)*dsin(gfl) 

do 78 i=l,n-l 
do 79 j=l,n 

plala2 = -p 1 +al *gt(j)-a2*dcos(gh)*gc(i) 
p2a2 = p2+a2*dsin(gh)*gc(i) 
ww = P 2a2**2+plala2**2 
wwm = P lala2**2- P 2a2**2 
if (m .EQ. 1) then 

parti = 2*dcos(2*gh)*plala2*p2a2**2-dsin(2*gh)*p2a2*wwm 
parti = partl/ww**2 

part2 = (-dsin(2*gh)*p2a2-dcos(2*gh)*plala2)/ww 

krl((k-l)*n+i,(l-l)*n+j)=(al/pi)*(partl+part2)*w(j) 

else 

if (m .EQ. 2) then 

parti = dcos(2*gh)*p2a2*wwm+2*dsin(2*gh)*plala2*p2a2**2 
parti = partl/ww**2 

kr2((k-l)*n+i,(l-l)*n+j)=(al/pi)*partl*w(j) 

else 

if (m .EQ. 3) then 

parti = -dcos(2*gh)*p2a2*wwm-2*dsin(2*gh)*plala2*p2a2**2 
parti = partl/ww**2 

part2 = (2*dsin(gh)**2*p2a2+dsin(2*gh)*plala2)/ww 

kr3((k-l)*n+i,(l-l)*n+j)=(al/pi)*(partl+part2)*w(j) 

else 

parti = 2*dcos(2*gh)*plala2*p2a2**2-dsin(2*gh)*p2a2*wwm 
parti = partl/ww**2 
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part2 = plala2/ww 

kr4((k-l)*n+i,(l-l)*n+j)=(al/pi)*(partl+part2)*w(j) 
end if 
end if 
end if 
79 continue 
78 continue 
end if 

311 continue 
310 continue 
320 continue 

********************* ******************** ********** 
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Appendix B - Schematic Form of Kernel Matrix [A] and Load- 
Vector {R} 
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Appendix C - FORTRAN Listing for Assemblage of Kernel Matrix 
[A] and Load Vector {R} 


The FORTRAN code for matrix [A] is 

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ** * * * * * * * * * * * * * * * * * 
* 

* Calculate the diagonal elements of matrix [A] 

* 

**************************************************** 

do 305 k=0,(2*ncr-l) 
do 50 j=l,n 
do 40 i=l,n-l 

a(k*n+i,k*n+j)=w(j)/(pi*(gt(j)-gc(i))) 

40 continue 

a((k+l)*n,k*n+j)=w(j) 

50 continue 
305 continue 

****************************************** ********** 
* 

* Combine krl,kr2,kr3 and kr4 to calculate matrix [A] 

* Input krl kr2 kr3 and kr4 to [A] 

* 

**************************************************** 

do 350 k=0,ncr-l 
do 351 l=0,ncr-l 
if (k .NE. 1) then 
do 352 i=l,n-l 
do 353 j=l,n 

a(2*k*n+i,2*l*n+j)=krl(k*n+i,l*n+j) 

a(2*k*n+i,2*l*n+n+j)=kr2(k*n+i,l*n+j) 

a(2*k*n+n-t-i,2*l*n+j)=kr3(k*n+i,l*n+j) 

a(2*k*n+n+i,2*l*n+n+j)=kr4(k*n+i,l*n+j) 

353 continue 
352 continue 
end if 

351 continue 
350 continue 

**************************************************** 
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The FORTRAN code for matrix [R] is 

****** ********************************************** 

* 

* Calculation of the right hand side of the systems equations R 

* 

**************************************************** 
do 370 i=l,ncr 

gsyiyi(i)=gs0xx*dsin(gf(i))**2+gs0yy*dcos(gf(i))**2 

gsyiyi(i)=gsyiyi(i)-gt0xy*dsin(2*gf(i)) 

gtxiyi(i)=-(gs0xx-gs0yy)/2*dsin(2*gf(i))+gt0xy*dcos(2*gf(i)) 

pixi(i)=-gsyiyi(i) 

qixi(i)=-gtxiyi(i) 

370 continue 
do 380 k=0,ncr-l 
do 381 i=l,n-l 

R(2*k*n+i)=-4*all*qixi(k+l) 

381 continue 
380 continue 

do 382 k=0,ncr-l 
do 383 i=l,n-l 

R(2*k*n+n+i)=-4*all*pixi(k+l) 

383 continue 

382 continue 

**************************************************** 
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Fig. 1. Geometry describing the interaction one notch with three microcracks. 
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Fig. 2. Mode- 1 normalized SIF, ki, for the inner crack tips versus the normalize 
d\2 — di 3 = <fn). Norihal far field stress, ayy, load case. 


0.5 


tip distance d (d = 


22 





Fig. 3 . Mode-II normalized SIF, k2, for the inner crack tips versus the normalize tip distance d 
(d = di2 = di3 = dn). Normal far field stress, <ryy, load case. 
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Fig. 4. Mode-I normalized SIF, kj, for the inner crack tips versus the normalize crack size a 3 (a 2 = 
a 4 = 0.1, constant distance between the center of each crack to the origin of the global 
coordinate system, |r,| = 1 ). Normal far field stress, ayy, load case. 
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Fig. 5. Mode- 1 1 normalized SIF, k 2 , for the inner crack tips versus the normalize crack size a 3 (a 2 - 
a 4 = 0.1, constant distance between the center of each crack to the origin of the global 
coordinate system, |r_,| = 1 ). Normal far field stress, <ryy, load case. 


25 



0 


0.2 


0.8 


1 


0.4 0.6 

Normalized crack size. a 4 


Fig. 6. Mode-I normalized SIF, kj, for the inner crack tips versus the normalize crack size a 4 (02 = 
03 = 0.1, constant distance between the center of each crack to the origin of the global 
coordinate system, |r,| = 1 ). Normal far field stress, <tyy, load case. 
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Fig. 7. Mode- 1 1 normalized SIF, k 2 , for the inner crack tips versus the normalize crack size a * (a 2 - 
a 3 = 0.1, constant distance between the center of each crack to the origin of the global 
coordinate system, |rj| = 1 ). Normal far field stress, cyy » load case. 
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Fig. 8. Mode- 1 normalized SIF, kj, for the inner crack tips versus the normalize tip distance d (d = 
d\2 = di3 = di4). Shear far field stress, <Jxy, load case. 
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Fig. 9. Mode-II normalized SIF, k 2 , for the inner crack tips versus the normalize tip distance d 
(d = dii = d X 3 = <fn). Shear far field stress, crxY, l° a< l case - 
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Fig. 10. Mode-I normalized SIF, ki, for the inner crack tips versus the normalize crack size <13 (aj = 
04 = 0.1, constant distance between the center of each crack to the origin of the global 
coordinate system, |rj | = 1 ). Shear far field stress, <txy, load case. 
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Fig. 11. Mode-II normalized SIF, k 2 , for the inner crack tips versus the normalize crack size a 3 (a 2 = 
o 4 = 0.1, constant distance between the center of each crack to the origin of the global 
coordinate system, |rj| = 1 ). Shear far field stress, <r X r, load case. 
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Fig. 12. Mode-I normalized SIF, kj, for the inner crack tips versus the normalize crack size (a 2 = 

<*3 = 0.1, constant distance between the center of each crack to the origin of the global 
coordinate system, Jr^- 1 = 1 ). Shear far field stress, ctxy , load case. 
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Fig 13 Mode- 1 1 normalized SIF, k 2 , for the inner crack tips versus the normalize crack s J ze ^ ( 3 
g ‘ « 3 = 0.1, constant distance between the center of each crack to the origin of the global 
coordinate system, |r,| = 1 ). Shear far field stress, a X y, load case. 


33 




Fig. 14. Geometry describing the interacting two notches with eight microcracks. 



Mode-1 normalized SIF, k 



Fig. 15. Mode-I normalized SIF, kj, for the inner crack tips versus normalized tip distance d (ai 
a 2 = 1, A = 0.1, a } = 0.1 for j = 3,. ..,10). 
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Fig. 16. Mode-II normalized SIF, k 2 , for the inner crack tips versus normalized tip distance d (oi 
a-i — 1, A = 0.1, a, = 0.1 for i = 3 101. 
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Fig. 18. Mode-II normalized SIF, k 2 , for the inner crack tips versus normalized offset A (aj = a 2 
d = 0.1, aj = 0.1 for j = 3,. ..,10). 
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Fig. 19. Mode- 1 normalized SIF, kj, for notches and cracks 4 and 9 tips and for the inner crack tips of 


cracks: 3, 5, 6, 7, 8, 10 (a, = a 2 = 1, A = 0.8, d = 0.1, a 3 = a s = a 6 = a 7 = a 8 = a 10 = 0.1). 


39 





Fig. 21. Mode-I normalized SIF, ki, for notches and cracks 4 and 9 tips and for the inner crack tips of 
cracks: 3, 4, 6, 7, 9, 10 (ai = a? = 1, A = 1.0, d = 0.1, a 3 = a 4 = a$ = aj = a 9 = ajo = 0.1). 
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Fig. 22 . Mode-II normalized SIF, k 2 , for notches and cracks 4 and 9 tips and for the inner crack tips of 
cracks: 3, 4, 6 , 7, 9, 10 (ai = aj = 1 , A = 1.0, d = 0.1, <* 3 = ola — = fl 7 = 09 = ajo = 0.1). 
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